/*
 * main.cpp
 *
 *  Created on: 26.02.2013
 *      Author: Roma Jam
 */

#include "fftw3.h"
#include "main.h"
#include "complex"

using namespace std;

#define SIZE    500

complex<double> Sin[SIZE],
                DftSin[SIZE],
                Sqr[SIZE],
                SpecSin[SIZE],
                DftSqr[SIZE],
                SpecSqr[SIZE];

fftw_plan dft_of_sin,
          dft_of_sqr;

double A = 1;
double freq = 1000;
double sampl = 0.01;

FILE *Txt, *W0;

int main(void){
    Txt = fopen("Sin.txt", "wt");

    for (int i = 0; i < SIZE; ++i) {
        Sin[i] = A * sin (2* (double)M_PI * freq * i);
        if( (i > 30) & (i < 45) ) Sqr[i] = 1;
        else Sqr[i] = 0;
    }

    dft_of_sin = fftw_plan_dft_1d((int)SIZE, (fftw_complex*)Sin, (fftw_complex*)DftSin, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(dft_of_sin);
    fftw_destroy_plan(dft_of_sin);

    dft_of_sqr = fftw_plan_dft_1d((int)SIZE, (fftw_complex*)Sqr, (fftw_complex*)DftSqr, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(dft_of_sqr);
    fftw_destroy_plan(dft_of_sqr);

    for (int i = 0;  i < SIZE; ++i) {
        DftSin[i] /= (double)sqrt(SIZE);
        DftSqr[i] /= (double)sqrt(SIZE);
    }

    RestoreSpectrum(DftSqr, SpecSqr, SIZE);
    RestoreSpectrum(DftSin, SpecSin, SIZE);

    printf("Save file");
    for (int i = 0; i < SIZE; ++i) {
        fprintf(Txt, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",
                real(Sin[i]), imag(Sin[i]),
                real(DftSin[i]), imag(DftSin[i]),
                real(Sqr[i]), imag(Sqr[i]),
                real(DftSqr[i]), imag(DftSqr[i]),
                real(SpecSqr[i]), imag(SpecSqr[i]),
                real(SpecSin[i]), imag(SpecSin[i])

        );
        printf(".");
    }

    printf("\nOk!\nFile had been saved\n");
    fclose(Txt);
}

void RestoreSpectrum(complex<double> *ArrIn, complex<double> *ArrOut, const int ASize){
    for (int i = 0; i < ASize; i++) {
        ArrOut[i] = (i < (ASize / 2))? ArrIn[(ASize / 2) + i] : ArrIn[i - (ASize / 2)];
    }
}
